Graded anharmonic crystals as genuine thermal diodes: Analytical description of 
rectification and negative differential thermal resistance 
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We address the heat flow study starting from microscopic models of matter: we develop an ap- 
proach and investigate some anharmonic graded mass crystals, with weak interparticle interactions. 
We calculate the thermal conductivity, and show the existence of rectification and negative dif- 
ferential thermal resistance. Our formalism allows us to understand the mechanism behind the 
phenomena, and shows that the properties of graded materials make them genuine thermal diodes. 
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Many works are devoted to the problem of understand- 
ing the heat flow starting from microscopic models of 
matter [ij , and most of them are carried out by means of 
computer simulations, sometimes with inconclusive re- 
sults. It creates a demand for analytical studies, but, 
since Debye, the microscopic models used to describe 
heat conduction are mainly given by systems of anhar- 
monic oscillators, which involve problems without precise 
solutions. Anyway, interesting properties have been dis- 
covered and their use proposed: e.g., the possibility to 
control the heat flow by using nano-devices such as ther- 
mal diodes, transistors, memories, etc 0-01 • The basic 
structure of these objects, the thermal diode, is a device 
in which heat flows preferably in one direction. There are 
analytical attempts to explain this phenomenon and/or 
design a diode by using simple methods [1, H, Q, but, 
again, most of the works are carried out by means of 
computer simulations [l-Q]. A recurrently used design 
of diodes is given by the sequential coupling of chains 
with different anharmonic potentials [2-01 ■ Although fre- 
quently investigated, it is criticized due to the difficulty to 
be constructed in practice 0. Recently, Chang et al. [7] 
built a diode in a experimental work by using a different 
procedure: graded materials, i.e., nanotubes externally 
and inhomogeneously mass-loaded with heavy molecules. 
Numerical computations [lo| also indicate rectification in 
a graded anharmonic system with abnormal conductivity. 

An important effect noticed in these studies is the neg- 
ative differential thermal resistance (NDTR) 0, [Hj, a 
phenomenon where the heat flux decreases as the ap- 
plied temperature gradient increases. NDTR is used to 
design a thermal diode with a big rectification factor; it 
is also crucial for the functioning of some models of ther- 
mal transistors and logic gates [l^. There are attempts 
to explain the origin of NDTR (in systems given by the 
coupling of different lattices), see e.g. the "phenomeno- 
logical approach" in ref . [13[ , but a general comprehensive 
understanding of the phenomenon is still lacking. 

Hence, considering these central subjects for the heat 
mechanism study, we address here the following issues: 
(i) the development of new methods of modeling the heat 
conduction problem in anharmonic systems; (ii) the an- 



alytical investigation of the graded mass system as a re- 
liable candidate for diode, different from that given by 
the coupling of different parts, whose rectification de- 
cays with the system size, and that is difficult to be con- 
structed in practise; (iii) the understanding of NDTR 
onset and related properties, in particular, in a nonlinear 
system that is not the coupling of different lattices. 

Here, we investigate graded anharmonic crystals with 
self-consistent reservoirs, details ahead, and show that 
graded materials are perfect candidates for diodes: their 
rectification does not decay with size (for certain mass 
distribution), they present NDTR, and may be con- 
structed in practice 0. We recall that rectification is 
absent in the classical harmonic version of this model 
14| . Our analytical formalism makes transparent the 
mechanism behind these phenomena. Rectification oc- 
curs because the total heat flow involves a sum of "local 
conductivities", see eqs.(|4l[5]), each one depending on the 
local temperature (for the anharmonic system, not for the 
harmonic one) and also on the masses of neighbor parti- 
cles. As we invert the system between two thermal baths, 
the distributions of masses and temperatures change in 
a different way, leading to a different heat flow - more 
comments ahead. For the NDTR onset, we have a com- 
petition between gradients of temperature and mass, see 
the denominator of the heat flow expression ([H [S]). For 
large gradients and anharmonicity, there is a change of 
the dominant term as we increase the temperature dif- 
ference, and so, NDTR appears. 

Let us introduce the model and our approach. We 
consider anharmonic crystals with stochastic reservoirs 
at each site. For simplicity, we take d = 1. We will work 
with the "self-consistent condition" , that means absence 
of heat flow between each inner reservoir and its site in 
the steady state, i.e., the inner reservoirs are not consid- 
ered as "real" thermal baths as those given by the reser- 
voirs at the boundaries: they describe only some resid- 
ual mechanism of phonon scattering not present in the 
Hamiltonian. The use of these hybrid models is recurrent 
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Precisely, we take N oscillators with Hamiltonian 



N 2 



where > 0, Jji — Jij, V is the anharmonic on-site 
potential: 'P{qj) — 9j/4; with time evolution 

1/2 

dqj ~ (pj /mj)dt, dpj = —{dH/dqj)dt—(jPjdt+'-fj dBj, 

(1) 

where Bj are independent Wiener processes; Q is the 
coupling between site j and its reservoir; and 7j = 
2QmjTj, where Tj is the temperature of the j-th bath. 
Here, we will study only nearest-neighbor interactions. 

The energy current inside the system is given by 
(J-j^), where (•) means the expectation with respect to 
the noise distribution, and 



Jj:3+i (qj ~ qj+i) 



Pj+i 
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(2) 



precisely, describes the heat flow from jth to (j + 
l)th site; details are found in previous works [lij . 

For ease of computation, we map our system on an- 
1, for all j. That is, we make the 
/rrTJqj, Pj — pj / Jm^ , and so. 



-1/2 



other with m. 



change of variables: Qj 



J, M and A are replaced by Jjk — {'nij)'~^l'^ J jk(mk) 
Xj ~ Xj/m^ Mj ~ Mj/mj. We will drop out the tilde 
notation in the unit mass system below, but we make the 
rescale later to come back to the general system. 

It is also useful to introduce the notation of the phase- 
space vector ip — {Q, P), with 2N coordinates. Then, the 
dynamics ([T]) becomes ip = —Aip — XV' {ip) + arj, where 
{A" + J) and a are 2N x 2N matrices 
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/ above is the unit N x N matrix; J is the N x N matrix 
for the interparticle interaction Jij] A4,T,T are diagonal 
N X N matrices: Mji = Mj6ji , Tji = QSji , Tji = 
TjSji. 1] are independent white-noises; V^f) is a 27V x 1 



matrix with V{(p)j = for j = I, . . . , N and V'{ip)i = 
dV{ipi-N)/dipi-N for i = iV+ 1, . . . , 2N. In what follows 
we use the index notation: i for index values in the set 
[N+l,N+2, . . . , 2A^], j for values in the set [1,2,..., N], 
and k for values in [1,2,..., 2A^]. 

In previous works [l^j we establish an integral repre- 
sentation for the correlation functions, and so, for the 
heat current, of systems with the stochastic dynamics 
considered here. It starts with a Gaussian measure, re- 
lated to the harmonic part of the interaction. Unfortu- 
nately, the analysis of the resulting formalism is still very 
intricate, in particular, for the case of hard anharmonic 
potentials. That is, it seems very difflcult to reach the 
anharmonic behavior starting from perturbations of the 
harmonic part of the system. Then, in other previous 
work 17f, we start an approximative scheme, that we 



conclude here, within this integral formalism in order to 
make it treatable. 

Let us describe our approach. Now, we first consider 
the equations of dynamics without the interparticle inter- 
action J, but with the anharmonic on-site potential. We 
do not know a strong solution for the decoupled anhar- 
monic problem , but we know the steady distribution: we 
follow Boltzmann, i.e., our system with J — involves 
only noninter acting particles, each one connected to a 
thermal bath, and so we have, in the notation Q, P-, 

N 

d^i,{Q,P) = exp( -Y,h\-'=''^ lTj)WdQjdPjlnoiui., 



H 



(J=0) 



iM.Ql + x,V{Q,) + \P] 



To turn on J, we use the Girsanov theorem, which re- 
lates the solution of the complete process ip (with J, the 
interparticle interaction) with the previous one <j) (with 
J — Q). Precisely, it states that , for ti,...,tk < t, 
{fri (ii) • • ■ t/frfc {tk)) = j (l>ri{h) ■ ■ ■ (j^r^ {tk)Z{t)diJ,, where 
(•) is the expectation for the complete process (p, d/x is the 
distribution associated to the expectations of (j) (the de- 
coupled process), and the "corrective" factor Z(t) is given 
by, after manipulations involving Ito calculus [l6|, ll7| , 



Zit) = exp {-ir^Mt)Jv(bj{t) + ir^MO)Jv^jiO)) exp i^j ds7~V«(s) Jij(/'j+jv(s)+ 

td3^,{s)jJ,jr^AUkis) ~ f^ds<j,,{s)jl^r^XV'{4>\{s) -\ ds<t,,{s)jl^r^J,,cty,{s)^ 



(3) 



We assume the boundary condition (/)(0) = 0, for simplic- u > N,v < N, i.e., J (j)u{t)(pv{t)Z(t)diJ,, etc. Writing 
ity. In the steady state, the heat flow is related to Z{t) = exp(— / W{(j)(s))ds, in a perturbative analysis, 
the expression limt^oo {'Puit)ipv{t) — ipu-N{t)(pv+N{t)) , we stay with terms such as J[4>u{t)(l)v{t)W{4>{s))]dsdfj,. 
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But we do not know the distribution d/i, that is very hard 
to calculate: for the nonlinear process we know only the 
steady distribution dn^, . Then, we introduce an approxi- 
mative scheme. 

First, to relate the fields (p{t) and (j){s), we use 
the Ito calculus which establishes that, for functions 
of ifim) = e-*«/((/.(0)), H - -^'^Vf + 
[A^(j) + XV' {(!))] ■ V , where V means the derivation in 
relation to (the index i, as well known, takes val- 
ues in [N + l,...,2iV]). The difference between the 
linear and nonlinear dynamics in the generator of the 
time evolution T-l above is in the term multiplying the 
gradient operator: precisely, instead of we have 
{A^ + XT"{(j))/(j))<j). Thus, to make easier the calculations, 
we replace (j)(t) by its average value. Moreover, in the ex- 
ponential relaxation of (p, we still replace V'{4))/(t) by its 
average value, more details ahead. All together means: 
(t){t) -> e~(*"^)'"(/)(s) e^(*"*^-^0(s), where A is given 
by A° with M replaced hy M = M + {XT" {(j)) / (/)) . We 
still have a problem: the computation of J 4>{s)4>{s)dfi is 
not possible, since we do not know the distribution dfi, 
as said before. Considering that we have an exponential 
convergence to the steady state, and so, the main terms 
involve s close to t, we propose to replace d/i by d/i,, the 
well known steady distribution. 

To summarize, our main approximations mean the re- 
placement of (l>{t) by {4>it)) and d/i by d/i*; after it, the 
expression for the heat flow will involve terms such as 

JdrJ dii, {e-^^4>) (e""-^</') W^(0), 



where the time dependence is carried only by exp(— r^), 
where r comes from t — s, and, as t ^ c», r € [0, oo]. 

In order to teste our approximative scheme, we first 
turn to the harmonic self-consistent chains, where rigor- 
ous results are known. For a system with particles with 
the same mass, and for the case of weak interparticle 
interactions, up to first order in J, we have 

lim {ifiu{t)'fiv{t)) = lim / (pu{t)(f>v{t)Z{t)dfi 

where r S [0, cxd]. After the r and (j) integrations, we 
get limt^^Juv {ipu{t)ipy{t)) = (J-™)2(2CM)-i(T„-r„). 
That is exactly the same value, considering the lower or- 
der in the interparticle interaction, of the rigorous com- 
putation [l^. This expression leads us to the correct 
thermal conductivity. Moreover, for the case of a chain 
with particles with alternate masses (two different val- 
ues), our scheme also works perfectly well: it gives, again, 
the same value of the rigorous computation. 

Let us, now, analyze our anharmonic crystal. Consid- 
ering first order in J, with the integration in r carried 
out after using a representation for e~'^-^ we get, for 
u> N,v <N, 



Cv)] X 



X {(M„7n ' + A/„7~') (0LiV0') JL + [An-iV {<f>u-NV'{cj,u-N)4>l) lu^ + A„ {<Pi^N'P'{cl^v)4>l) 7."'] JL] , (4) 

I 



where Mu = Mu^n, Duv = [Mu - Mvf + {MuCv + 
MvCu){Cu + Cf )■ For u> N, {(1)1) = r„; but the compu- 
tation of {(f>1), V < N, is not easy (note that d^* is a sin- 
gle variable distribution, and so, {4>u4'™) — (^m) (</>™))- 
Let us assume some regime before any approximation: 
we consider a high anharmonic system, i.e., A large 
and M small. Thus, we take {4>l) = Ic^tV^ j^^J"^ , 
{(j)^) = 4c2Ti,/A„. If M = 0, we would have ci ~ 
r(3/4)/r(l/4) ~ 1/3, C2 ~ r(5/4)/r(l/4) = l/4. To de- 
termine the values of ci and C2 , we turn to the expression 
of the heat current = Juv{{(l>u4'v) ~ {(t>u-N(t>v+N)) f^, 
with u — N — j, V = j + 1, take all sites at the same 
temperature T and find the values such that J^j^ = 0. 
We obtain C2 — 1/4 and ci — 1/2. Then, we per- 



form the further computations. For high anharmonicity 
and very small temperatures, for the dominant term in 
J-^j^ = J^jj+i, we obtain, after the rescaling back to the 
system with general mass values, i.e., Xj — )■ Xj/m3, etc. 



-1 



(5) 



1/2 



where we take, after the rescale, uniform potentials and 
couplings: Xj = X, etc. From J^j.j+i above and the self- 
consistent condition T = T\^2 — ^za = ■ ■ ■ = Tn-i,n, 
which establishes that the heat current comes from the 
first reservoir, passes through the chain and goes out by 
the last reservoir, we determine the temperature profile. 
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We have 

Tim2T^^^ + miT^^^)/C = Ti - T2 

T{mNTl/^^ + mN^iTl/^)/C = Tn-i - Tn, 
where C — J^/\^^^(^. We sum ah the equations to obtain 
J- {(7712^1 + (mi + m3)T^^^ + ... 

+ {mN-2 + mN)Tl/\ + TOAr_iry^| /C = Ti - Tat, 

that gives us, from = JC{Ti — Tn)/{N — 1), an ex- 
pression for the thermal conductivity /C. The system of 
equations above may be rewritten as 

For the case of particles with the same mass, the equa- 
tions become T^^^ - T^^^ = . . . = tI/\ - T^^^ that 
leads to a linear profile for T^/^, i.e., T^^^ = T^^ + [{k - 

1)/{N - 1)]{tI/^ - T//^). For a general mass distribu- 
tion, the problem is more complicated: let us examine 
it in the case of a small temperature gradient. We write 
Tl = T + aie and Tpf — T + aNe; T, ai, un and e given (e 
small). Then, Tj is a function of e, with values between 
Tl and Tat: Tj = 7+0^6 + ©(e^). Let us analyze only the 
first order in e. From the equations for the self-consistent 
condition, we get the solution aj = ai + (ai — aAr)5j/iSjv, 
Sj — mi + 2m2 + . . . + 2toj_i + rrij. Hence, turning to 
the thermal conductivity formula, after algebraic manip- 
ulations, we obtain 

j_ _ J_ _ £(ai - gjv) r 2 _ 2] 
IC IC' ^ CiN - 1)2T^/^Sn ^"^^ ""'J' 

where /C' is the conductivity for the system with inverted 
boundary baths. And so, there is rectification even for a 
small gradient of temperature. And more, if the graded 
mass grows with iV^, i.e., rrij = p ■ mi, then the differ- 
ence above does not decay with N . By taking T^q > Ti 
(i.e., ajv > 0,1) and mi^ > mi, we see that the thermal 
conductivity is bigger when heat flows from the large to 
the small mass, as experimentally observed in a graded 
system 0. 

We stress here that the dependence on temperature 
for the local anharmonic conductivity comes from the 
dynamics: {4>{t)) ~ e^~*^^(/)(0), where A depends on T 
for the anharmonic (not for the harmonic) case. The 
combination of particle masses and temperatures, and 
the difference as we invert the chain, lead to rectification. 

Now, we consider the investigation of NDTR. We turn 
to eq.(j4]), that is directly related to the heat flow, and is 
valid for weak interparticle interaction J in any regime: 
low and high anharmonicity, temperature, etc. All the 



terms include D in the denominator, except the first 
one that, however, may be manipulated and absorbed 
by the other terms. Hence, will have Djjj^i in 

the denominator (see e.g. the first equality in eq.®, 
the expression for high anharmonic regime), where, we 
recall, Dj^j+i = {Mj - Mj+i f + 2(^^{Mj + Mj+i), 
A4j = Mj + (Aj0|), expression determined for a sys- 
tem with unit masses. For high anharmonicity we have 
(Xj4>'^') ^ T^^'^ / y^J"^] and for very low anharmonicity, 
(Xjcj)^) ^ XjTj/Mj. Rescaling to get the expression for a 
system with different values for the particle masses, and 
considering high anharmonicity, just to fix the expression 
for the temperature behavior (but the analysis below, ad- 
justing the power of T, follows anywhere), we have, for 
the first term in Dj^j^i 

{Mj - Mj+if = [cjmj+i - Cj+imjf /{mjmj+if, 

Cj = [AI + X^/^T^^]. The second term in D always in- 
creases with T, and is subdominant for ^ small: precisely, 
for < X'-^'^ AmjT^^^ /mjmj+i; note however that it 
shall dominate for very small T, as assumed in the sec- 
ond part of eq. (O . Let analyze the first term, considering 
a graded mass chain. For m^+i > mj and T^+i > Tj, i.e., 
gradient of mass and temperature at the same direction, 
if ATj = Tj^i — Tj <^ Amj = mj^i — mj, then Cjmj^i > 
Cj+imj-; and if ATj 3> Am-,-, then Cj-mj+i < c^+im^-. Re- 
call that ATj increases as we increase Tj, and, of course, 
ATj depends also on Ti — Tn , the "total gradient" : ATj 
will be very small if |Ti — Tn\ is very small. Hence, 
starting from a very low total temperature gradient AT, 
as we increase AT, then {Mj — Mj+i)"^ first decreases, 
but after same point it becomes an increasing function. 
That is, 1/D first increases and, in sequel, decreases with 

AT. As we have T — J-j ^ D~^ATj, and D changes as 
~ 1/2 . ~ . 

cATj , with c depending on A, Am^, if A and Ato^ are 

1/2 

not very small, then cAT^ dominates ATj (Tj < 1), 
and the heat current first increases and then decreases 
with AT. In other words, we have NDTR. 

To conclude, we stress that diodes of graded materials 
sound to be experimentally reliable 0, and ubiquitous 
structures: our results follow for many other anharmonic 
potentials as indicated by the formalism derivation. 
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